Continuous thermal melting of a two-dimensional Abrikosov vortex solid 
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We examine the question of thermal melting of the triangular Abrikosov vortex solid in two- 
dimensional superconductors or neutral superfluids. We introduce a model, which combines lowest 
Landau level (LLL) projection with the magnetic Wannier basis to represent degenerate eigenstates 
in the LLL. Solving the model numerically via large-scale Monte Carlo simulations, we find clear 
evidence for a continuous melting transition, in perfect agreement with the Kosterlitz-Thouless- 
Halperin-Nelson- Young theory and with recent experiments. 



The problem of melting of two-dimensional (2D) solids 
has been around for several decades. As was first pointed 
out by Kosterlitz and Thouless a 2D solid to liq- 
uid transition can generally be expected to be very dif- 
ferent from its three-dimensional (3D) counterpart. 3D 
solids melt via a first-order transition, where the magni- 
tude of the long-range crystal order, measured e.g. by the 
strength of the Bragg peaks in the structure factor, drops 
to zero discontinuously from a finite value just below the 
melting temperature T m . In 2D, long-range crystal or- 
der, which breaks continuous symmetry of spatial trans- 
lations, is impossible at any nonzero temperature @]. 2D 
solids are thus characterized by power-law decay of crys- 
talline correlations at low temperatures. This, however, 
is enough to give rise to a nonzero shear modulus and the 
low temperature phase is thus a true solid. The order in 
this case is topological, in the sense that the solid phase is 
characterized by the absence of free topological defects, 
i.e. dislocations, which are bound into pairs with Burg- 
ers vectors equal in magnitude and opposite in direction. 
Kosterlitz and Thouless proposed [l[ that melting in 2D 
can happen continuously via the unbinding of dislocation 
pairs. Such a melting transition is then closely analogous 
to the well-known Kostcrlitz-Thouless transition in 2D 
superfluids. Halperin, Nelson and Young Q developed 
the idea of Kosterlitz and Thouless into a detailed the- 
ory of dislocation-mediated 2D melting transition, which 
is now frequently referred to as KTHNY theory 

Experimental confirmation of KTHNY theory has 
proven to be somewhat difficult to obtain [4(. However, 
by now there exist reports in the literature of appar- 
ently continuous melting transitions of triangular solids, 
which agree very well with KTHNY theory predictions 
[H. In particular, very recently direct STM imaging of 
dislocation-mediated melting of a triangular vortex solid 
in a thin-film superconductor has reported a continuous 
transition |(|. In contrast, there exists no direct theoret- 
ical evidence of a continuous finite-temperature melting 
transition in any microscopic model of 2D vortex solids 

In this paper we address the problem of the melt- 
ing of Abrikosov vortex lattices in 2D superconductors 
and neutral superfluids by introducing a model which 
combines lowest Landau level (LLL) projection with the 



magnetic Wannier basis to represent degenerate eigen- 
states in the LLL. Solving the model using state of the 
art Monte Carlo simulations, we obtain clear evidence 
for a continuous melting transition, in perfect agreement 
with KTHNY theory. 

We start from the standard Ginzburg-Landau (GL) ex- 
pression for the energy functional of a superconductor: 
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The superconductor is placed in a perpendicular mag- 
netic field B = V x A. We will assume that we are deal- 
ing with a strongly type-II superconductor, as is the case 
in the experiment of Ref. Q. In this case fluctuations 
of the electromagnetic field can be neglected and our re- 
sults will then be applicable to neutral 2D superfluids 
in an artificial perpendicular magnetic field as well. The 
quadratic part of Eq. ([1]) is diagonalized by expanding the 
complex order parameter ^(r) in terms of Landau level 
eigenstates. If the magnetic field is close to the upper 
critical field H C 2, or, equivalently, the superconducting 
coherence length is of the order of the magnetic length 
I = hc/eB, the lowest Landau level (LLL) approxima- 
tion can be used [1, Q, when only the contribution of 
the LLL to the eigenstate expansion of the order param- 
eter is retained. We will adopt the LLL approximation 
henceforth. The order parameter is then written as: 



i0m(r). 



(2) 



Here m is as yet unspecified LLL eigenstate label, m (r) 
is the corresponding cigenfunction and c m is the complex 
amplitude, corresponding to the LLL eigenstate m (r). 
A crucial ingredient of our work is a judicious choice 
of the LLL basis in Eq.([2]), which should be chosen in 
such a way as to lead to the simplest representation of 
the low energy states of the system. As will be clear 
from the discussion below, the best LLL eigenstate basis 
for our problem is the basis of magnetic Wannier func- 
tions, first introduced in Ref. (lp| . and employed e.g. in 
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111 1 1 21 ] . We will only briefly mention the properties of 
magnetic Wannier functions that will be important for 
our present problem. Readers interested in a more de- 
tailed discussion should consult Rcfs. Magnetic 
Wannier states are defined on the sites of a triangular 
lattice with basis vectors ai = ax, a 2 = a(x + y/3y)/2, 
so that r m = miSLi +TO2&2 with integer mi : 2- The lattice 
constant a 2 = 47r£ 2 /y/3 is chosen in such a way that the 
unit cell of the lattice contains exactly one magnetic flux 
quantum |ai x a2| = 27r£ 2 . The explicit form of magnetic 
Wannier functions 4> m (r) will be unimportant for us, we 
will only mention that these states are normalizable but 
have a 1/r 2 decay at long distances, which is the fastest 
decay compatible with the LLL projection. The utility 
of this set of LLL eigenstates for our problem becomes 
apparent when one notices that, unlike arbitrary spatial 
translations in the presence of a perpendicular magnetic 
field, translations by lattice vectors of the above lattice of 
magnetic Wannier states commute with each other. This 
allows us to define lattice momentum k, belonging to the 
first Brillouine zone of the triangular lattice and the cor- 
responding Bloch states vE'k(r) = Yl m 0me lk rm , 
where Nj, is the total number of ma gne tic flux quanta 
piercing the sample. As was shown in [12J , the wavefunc- 
tions ^(r) are nothing but the magnetic Bloch states, 
first introduced by Eilenberger [l3| . These states form a 
complete orthonormal set of states in the LLL and cor- 
respond simply to triangular Abrikosov vortex lattices, 
with vortex cores located at: 



r m + (ai + a 2 )/2 + ^ 2 z x k, 



(3) 



i.e. the Bloch momentum k labels the center of mass 
positions of different Abrikosov lattice states. As shown 
in Ref. [ijj], the GL functional Eq. (Q} takes a particu- 
larly simple form, when written in the magnetic Wannier 
basis. Namely, the quartic term of the GL functional 
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where I mim2m3m4 = / d 2 rcf>* mi (r)0 m2 (r)^ m3 (r)0 mj (r), 
turns out to possess a low-energy symmetry correspond- 
ing to center of mass conservation, i.e. imim 2 m 3 m 4 ~ 
<5mi+m 2 ,m 3 +m 4 - This symmetry is a reflection of the 
translational symmetry of the 2D plane, leading to the 
degeneracy of all the Abrikosov vortex lattice solutions, 
and of the LLL projection, which makes the set of 
Abrikosov vortex states a complete set. Taking into ac- 
count that the matrix elements / m im 2 m 3 m 4 are short- 
range, namely have a 1/r 6 decay at large distances, we 
then arrive at the following LLL representation of the GL 
functional: 



0, 



(4) 



where P labels all possible smallest 4-site plaquettes of 
the triangular lattice (see Fig. 1) and we have taken 
c m ~ e m , while |c m | 2 is assumed fixed (note that this 



is not the same as neglecting fluctuations of the ampli- 
tude of vE'(r), which would lead to Landau level mixing). 
K can be easily related to the parameters of the origi- 
nal GL functional Eq.([T]), but will be left here as a phe- 
nomenological parameter. Its physical meaning, as will 
become clear shortly, is the shear modulus of the vor- 
tex lattice. Eq.Q represents the shortest-range phase- 
dependent center of mass conserving quartic term on the 
triangular lattice. The quadratic term in the GL func- 
tional and the phase-independent quartic terms simply 
determine the magnitude of |c m | 2 and consequently of 
the parameter K. 

Let us now demonstrate that Eq.Q indeed properly 
represents elasticity theory of the Abrikosov vortex solid. 
It is easy to show by direct substitution that the set of 
minimum energy states of Eq.Q corresponds to all pos- 
sible states with uniform gradients of the phase 9 m along 
the basis directions of the triangular lattice [lij ]. Sub- 
stituting the corresponding amplitudes c m ~ e 10 ™ 1 into 
Eq.@ one obtains precisely the set of magnetic Bloch 
states <3/k, where each momentum k corresponds to a 
given value of the phase gradient. It is then clear from 
Eq. ([3]) that in the long- wavelength elasticity theory of the 
vortex lattice, we can identify local vortex displacements 
with gradients of the phase 0: 



u = tzx ve. 



(5) 



Defining the strain tensor Uij = {diUj + djUi)/2 in the 
standard way and taking into account the incompressibil- 
ity of the vortex lattice Uu = (summation over repeated 
indices is implicit), which immediately follows from 
one obtains the following expression for the clastic en- 
ergy: 



E 



d 2 ruj 
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where (i is the shear modulus of the vortex lattice. To 
connect Eqs.Q and (|6|), we expand the cosine in Eq.(|4]) 
to quadratic order in the phase, discard the ground state 
energy term, and take the continuum limit: 
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9K a 4 
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(7) 

which gives K = 2/i£ 2 /37r and defines the physical mean- 
ing of K. The lack of the (V6*) 2 term in (|7|) is a conse- 
quence of the center of mass conservation symmetry of 
Eq.©. Eq.(j4|) can thus be thought of as a lattice rcgular- 
ization of the continuum elasticity theory of the incom- 
pressible Abrikosov vortex solid Eq. (|6]) . The representa- 
tion of the continuum elasticity theory of the Abrikosov 
vortex lattice in terms of the phase Laplacian Eq. ([6]) has 
in fact been obtained before by Moore [HI]. It has not, 
however, been realized that a lattice rcgularization of ([7]l. 
Eq.(HJ), can be used to study the melting transition. 
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Before turning to numerical investigation of the phase 
diagram of the lattice model Eq. Q, let us calculate 
some of the basic properties of the vortex solid phase 
using the continuum theory ([7]). In particular, corre- 
lations in the vortex solid phase can be quantified by 
two correlations functions: the phase correlation func- 
tion G p h(r) = (^e ie ^e~ l6 ^) and the vortex position cor- 
relation function G v (r) = ( e *G-[u(r)-u(o)]^ where q ig 

a reciprocal lattice vector. Straightforward calculation 
gives: 



G v (r) ~ 1/r" 



Vg 



\G\ 2 T 



(8) 



while G p h(r) vanishes in the thermodynamic limit for 
any nonzero r and at any nonzero temperature. The 
exponent rjc is in perfect agreement with the KTHNY 
results Q . The short-range nature of the phase correla- 
tion function G„h in the vortex solid has been emphasized 
by Tesanovic jig , and our results are in agreement with 
Ref. 0. 

As discussed above, we expect the power-law vortex 
solid order Eq. © to exist up to a melting temperature 
T rn , at which free dislocations will appear and make the 
vortex positional correlations short range. To calculate 
T m we will employ the standard argument [l| , comparing 
energetic and entropic contribution of a single dislocation 
to the free energy of the system. We assume the presence 
of an isolated dislocation with Burgers vector b = ax: 



Vw x ■ d£ = a, 



(9) 



where the contour in the above integral encloses the dislo- 
cation core. Solution for the displacement field, minimiz- 
ing the elastic energy of an incompressible solid Eq. © 
under the constraint ©, is given by u x = ac^jl-n + 
(a/4ir) sin(2(ft), u v = — (a/4n) cos(20), where 4> is the 
polar angle [17]]. Substituting this solution into Eq.©, 
one obtains the following result for the dislocation energy 
Ed = (^a 2 /2tt) ]n(L/a), where L is the size of the system. 
Comparing this with the entropic contribution of the dis- 
location to the free ener gy — T\n(L/a) 2 , one obtains the 
melting temperature 0, 17 1 T m = fia 2 /4ir = \/3t:K/2, 
which gives: 



K/T m = 2/V3n. 



(10) 



This can be compared with our numerical results, if K is 
replaced by the actual value of the shear modulus at T m 
(see Eq.(TTT]) below). 

According to KTHNY theory Q, the dislocation- 
mediated melting of a triangular solid results in an in- 
termediate hexatic phase, between the solid and the 
isotropic liquid, which has power-law orientational order 
of the solid, but no positional order. We do not expect 
to observe this phase in our numerics, which we will de- 
scribe shortly. The reason is that our choice of the LLL 
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FIG. 1: (Color online) The internal energy, (H), of Eq. (gj, 
using Monte Carlo simulations for two different system sizes. 
The inset shows the labeling of sites on the three different 
plaquette orientations in Eq. Q. 



basis and periodic boundary conditions that this choice 
implies, fixes the orientation of the vortex solid, thus pre- 
cluding spontaneous orientational order. This is true for 
the simulations of Ref. [1] , which used the Landau gauge 
orbital basis, as well. We then expect to observe a single 
solid to liquid transition. The main question we would 
like to resolve is whether the transition is driven by the 
unbinding of dislocations and described by KTHNY the- 
ory, or it is first order as in the simulations of Ref. Q 
and thus possibly unrelated to unbinding of topological 
defects entirely 

We now turn to numerical Monte Carlo simulations of 
our lattice model Eq. ((4]). Taking the phases 9 m to be 
continuous variables between and 2-7T, we use a modified 
Metropolis algorithm that attempts an update on each 
phase variable selected randomly from a uniform distri- 
bution between and A# max . An important observation 
we make is that the choice of A6* max critically affects 
the ergodicity of the update. In conventional simula- 
tions of XY-type models, A# max is varied as a function of 
temperature in order to maximize simulation acceptance 
rates. We find that in our model this leads to problems 
with freezing into nearby metastable states. In order to 
combat this, we systematically explored the simulation 
dynamics as a function of a temperature-independent 
A# max . We find that for insufficiently large A6* max the 
simulation freezes into a metastable state of higher en- 
ergy in an intermediate temperature regime, however, 
with sufficiently large A# max , the simulation is ergodic, 
and always finds the correct minimum free-energy state 
[l9j | . Figure [T] shows the simulation energy in the vicinity 
of a finite-temperature phase transition (T m « 1.3-KT, dis- 
cussed below) for two different system sizes. No evidence 
of a first-order discontinuity (latent heat) is apparent. 

Figure [2] shows the main result of our simulations, 
which is the finite-size scaling behavior of the shear mod- 
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FIG. 2: (Color online) The finite-size scaling behavior of the 
shear modulus through the finite-temperature transition. The 
lower-left inset is data for linear lattice size L = 16 over a 
larger temperature range. The solid black line is the equation 
K = (2/V~3n)K/T, which from Eq. (TO} gives a finite-size 
estimate for T m (L). The upper right inset is the finite-size 
scaling [T^ of T m (L), as discussed in text, which gives an 
estimate for T m (oo) (the star) of 1.32 K. 



ulus JC, defined as the second derivative of the free energy 
with respect to the uniform shear angle tp s = \/37r£ 2 V 2 9: 




where the prefactor in the expression for ip s is chosen so 
that K(T = 0) = K . As is apparent from Fig. [2] we ob- 
serve a single continuous phase transition from the vortex 
solid state with a nonzero shear modulus, to the high- 
temperature liquid phase. Comparing the value of the 
ratio IC(T m )/T m with the value in Eq. (fTUJ), we see per- 
fect agreement with the KTHNY theory. We can make 
a highly-accurate quantitative assessment by compar- 
ing the finite-size scaling behavior of T rn obtained using 
Eq. (flUI) to the expected scaling form T m (L) — T m (oo) oc 
T m (oo) ln(L)- x / Q (Fig. H inset), where a = 0.36963, 
from which we extract T m (po)/K sa 1.32. This scaling 
form follows from the KTHNY expression for the tem- 
perature dependence of the positional order correlation 
length £ ~ exp[C/(T - T m ) a ] |. 

In conclusion, we have introduced a microscopic model 
for the melting of a 2D vortex solid using the magnetic 
Wannier basis to represent degenerate eigcnstates in the 
LLL. Solving the model using large-scale Monte Carlo 
simulations shows the finite-temperature melting tran- 
sition is continuous, in agreement with KTHNY theory 
and recent experiments. We note that previous work on 
related model have observed a first-order melting transi- 
tion Q. We believe that the most likely reason for this 
difference may be our neglect of the fluctuations of |c m | 2 , 
only retaining fluctuations of the phase 6 m (note again 
that this does not mean we are neglecting fluctuations of 
the amplitude of the order parameter Vf^r)). The |c m | 2 



fluctuations, while not soft, could conceivably lead to a 
weak fluctuation-induced first order transition. An indi- 
rect confirmation of this scenario is the fact that the melt- 
ing temperature in Ref. Q is observed to be very close 
to the KTHNY melting temperature. Therefore, our re- 
sult clearly demonstrates that dislocation unbinding is 
the mechanism behind the vortex lattice melting tran- 
sition in 2D, even though the transition may be driven 
weakly first order in specific cases. 
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